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Abstract — Many biological oscillators are arranged in net- 
works composed of many inter-coupled cellular oscillators. 
However, results are still lacking on the collective oscillation 
period of inter-coupled gene regulatory oscillators, which, as has 
been reported, may be different from the oscillation period of an 
autonomous cellular oscillator. Based on the Goodwin oscillator, 
we analyze the collective oscillation pattern of coupled cellular 
oscillator networks. First we give a condition under which 
the oscillator network exhibits oscillatory and synchronized 
behavior, then we estimate the collective oscillation period based 
on a multivariable harmonic balance technique. Analytical 
results are derived in terms of biochemical parameters, thus 
giving insight into the basic mechanism of biological oscillation 
and providing guidance in synthetic biology design. Simulation 
results are given to confirm the theoretical predictions. 

I. Introduction 

Diverse biological rhythms are generated by multiple 
cellular oscillators that somehow manage to operate syn- 
chronously. In systems ranging from circadian rhythms to 
segmentation clocks, it remains a challenge to understand 
how collective oscillation patterns (e.g., oscillation period, 
amplitude) arise from autonomous cellular oscillators. As 
has been reported in the literature, there can be significant 
differences between collective oscillation patterns and cell 
autonomous oscillation patterns, the difference are embodied 
not only in the oscillation amplitude [1], but also in the 
oscillation period [2], [3]. 

The famous Goodwin oscillator provides a perfect model 
to study the mechanism of how collective oscillation pattern 
arises from autonomous cellular oscillators. The Goodwin 
oscillator was proposed in 1965 to model the oscillatory 
behavior in enzymatic control processes [4], [5]. It is a 
minimal model that describes the oscillatory negative feed- 
back regulation of a translated protein which inhibits its own 
transcription. Because the Goodwin oscillator can capture 
the essential characteristics in biochemical oscillators, it has 
been extensively used to model biological oscillators such 
as ultradian clocks of vertebrate embryos [6] and circadian 
clocks of neurospora [7], drosophila [8] as well as mammals 
[9]. 
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Another advantage of the Goodwin oscillator is that it 
allows for an analytical understanding of the basic dynamical 
mechanisms, whereas other biophysically substantiated mod- 
els of cellular oscillators usually exhibit high complexity and 
large number of variables, which hamper the mathematical 
treatment, and subsequently obscure the underlying mecha- 
nisms. The Goodwin oscillator allows one to gain insights 
into the mechanisms of biochemical rhythms. For example, 
the oscillation conditions of a single Goodwin oscillator 
were obtained in [10], [11], [12], [13]. The synchroniza- 
tion conditions for coupled Goodwin oscillator networks 
were reported in [14], [15], [16]. Based on a multivari- 
able harmonic balance technique, the oscillation patterns of 
a single Goodwin oscillator were obtained in [17], [18]. 
This is an important step toward understanding the period 
determination in biochemical oscillators. However, given 
that biological rhythms are generated by multiple cellular 
oscillators coupled through intercellular signaling, it remains 
a challenge to determine the periods in biological rhythms, 
which ranges from seconds in cardiac cell contraction to 
years in reproduction. Recently, using the phenomenological 
phase model, the authors in [19] proved that if the inter- 
cellular coupling is weak, the collective period is identical 
to the autonomous period. However, since the phase model 
contains no direct biological mechanism of the cellular clock, 
it can potentially weaken the model's reliability in checking 
scientific hypotheses. 

This paper derives analytical results for the collective os- 
cillation period of multi-cellular networks based on coupled 
Goodwin oscillators. Specifically, we study the collective 
period of a network of Goodwin oscillators connected by 
diffusive coupling. The basic idea is to use a multivariable 
harmonic balance technique [20]. In fact, due to the multi- 
cellular structure, the solution to harmonic balance equations 
in the multivariable harmonic balance technique is very 
difficult to obtain. Here we are interested in the collective 
period, so we can circumvent the problem by restricting our 
attention to solutions corresponding to synchronized oscilla- 
tions in the oscillator network. To this end, we also give 
an oscillation/synchronization condition for the Goodwin 
oscillator network. 

II. Model description and transformation 

The Goodwin oscillator describes the dynamics of an os- 
cillatory negative feedback regulation loop. As shown in Fig. 
[H mRNA Xi produces a protein X2 which, in turn, activates 
a transcriptional inhibitor X3. The inhibitor X3 inhibits the 
production of mRNA Xi, which closes a negative feedback 
loop. The kinetic dynamics of the Goodwin oscillator are 



given by [4], [10], [21]: 
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Here [Xi], [X2], and [X3] are concentrations of mRNA Xi, 
protein X2, and inhibitor X3, respectively; vq, vi, and V2 are 
the rates of transcription, translation, and catalysis; ki, fc2, 
and k^ are rate constants for degradation of each component; 
I /Km is the binding constant of end product to transcription 
factor; and p is Hill coefficient, which describes cooperativity 
of end product repression. 
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Fig. 1. Schematic diagram of a Goodwin oscillator. The mRNA Xi 
produces a protein X2 which, in turn, activates a transcriptional inhibitor 
X3. The inhibitor X3 inhibits the production of Xi, wliich closes a negative 
feedback loop. 



There are other interpretations of the above Goodwin 
oscillator [4], for example, one may regard X2 as an enzyme 
precursor that after primary synthesis on mRNA templates, 
Xi, passes through a pool of inactive molecules before being 
transformed into mature, active enzyme, X3. One may also 
take Xi to be an enzyme population whose rate of synthesis 
is regulated by feedback control at the polysome level via 
a metabolite X3. In this case, X2 is then an intermediate in 
the biosynthetic sequence leading to X3. 

The Goodwin oscillator ([T) can be transformed into 
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via the introduction of dimensionless variables 
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X2 = -^ , x^ = -—, t = -T 

In (|2]l, 61, &2, and 63 are positive parameters given by 

hi = ki<;, i = 1,2,3 

Suppose the oscillator network is composed of N oscilla- 
tors, and they are connected by diffusive coupUng, then the 



dynamics of the network are given by 
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where i = 1, 2, . . . , iV denotes the index of the ith oscillator, 
and Qij > denotes the coupling strength between oscillator 
i and oscillator j. If Uij = 0, then there is no interaction 
between oscillator i and oscillator j. 

Assumption 1: We assume that the interaction is bidi- 
rectional, i.e., ai,j = ttj^i. We also assume that 
the interaction topology is connected, i.e., there is a 
multi-hop path (i.e., a sequence with nonzero values 

j) from each node i to 



''2,rni 5 ^mi ,7712 : 
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every other node j. 

Remark 1: Since protein is usually diffusible and involved 
in intercellular signaling [22], we model the intercellular 
coupling by the diffusion of X2, which usually represents 
the protein product in the Goodwin oscillator model. 

For convenience in analysis, we can write (|3]l into the 
following matrix form: 
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Next, based on a multivariable harmonic balance tech- 
nique, we will analyze the collective period of the Goodwin 
oscillator network in (|4). 

III. Oscillation/synchronization condition 

To study the collective period, we need to guarantee that 
the component elements in Xi {i = {1,2,3}) in (|4|i first 
oscillate, and then further more, oscillate in synchrony. In 
this paper, we consider the Y-oscillation, which is defined as 
follows [20]: 

Definition 1: A system x = f{x) with x{t) G 7?,'" is said 
to be Y-oscillatory if each solution is bounded and there 



exists a state variable Xi such that lim Xjjt) < lim Xi{t) 

for aknost all initial states a;(0). 

To prove that (|4]i is Y-oscillatory, we introduce Lemma [T] 

Lemma 1: [23] System ^ is Y-oscillatory if all condi- 
tions (a), (b), and (c) hold; 

(a) It only has isolated equilibria X*; 

(b) The positive semiorbit {X{t) = 
[Xi(t),X2{t),X3{t)]\t > and f e domX(«)} is 
bounded; 

(c) The Jacobian matrix evaluated at the equilibrium point 
X*, i.e., J, has at least one eigenvalue with positive 
real part. 

Proof: The Lemma can be obtained by combining 
Theorem 1 in [23] and the discussion below its proof which 
shows that the hyperbolicity condition can be relaxed. ■ 

In the following, we will prove that ^ satisfies all the 
conditions in Lemma [T] and hence is Y-oscillatory. We will 
also give a synchronization condition for the oscillation. In 
this manuscript, we define synchronization as follows; 

Definition 2: System ^ is synchronized if the elements in 
X3 satisfy lim \x3 i{t)—X3_j{t)\ = for any 1 < i, j < N. 

Remark 2: Only xs^i [i ~ 1,2, ... ,N) is used in the 
definition of synchronization. This is because according 
to the modeling assumption, x^^i is corresponding to the 
concentration of inhibitor or enzyme, which can be regarded 
as the output of an Goodwin oscillator Note that under 
this definition, when the system is synchronized, xi ^ {i — 
1,2, ... ,N) may be identical or non-identical. The same 
situation holds for X2,i {i = 1,2,..., N). 

Theorem 1: Q has oscillatory solutions if it satisfies the 
following inequality 
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(61 + &2 + ^3)(&1&2 + ^2&3 + ^1^3) - &1&2&3 



>1 (7) 



where xq is the unique positive solution to 1/(1 + Xg) = 
h\h2hi,XQ. Furthermore, the oscillations in all oscillators are 
synchronized if the algebraic connectivity q (which is defined 
as the second smallest eigenvalue of matrix A in (|6]l) of 
coupling topology satisfies 



7 f V^ 
Q > ^oi + ,, , ; 7 = max <^ -; — 



(8) 



Proof: From Lemma [T] we know to guarantee a Y- 
oscillatory solution, we need to prove the conditions (a), (b), 
and (c), are satisfied. The proof is decomposed into three 
steps. 

Step I — Satisfaction of condition (a): Using the mono- 
tonic property of /(•), we can prove that (O has only 
one equilibrium point X^ = [ xq xq ... xq ] with 
xq > determined by /(a;o) = 616263X0. The derivation is 
as follows; 

In the equilibrium point, we have 



dxi 

dt 



0,* = {1,2,3} 



So after some simple algebra, (|4]i is reduced to 

/(X*) - 616263X3* - 6163AX3* (9) 

where /(X3) and A are defined in (|5]l and (|6]l, respectively. 
To prove the uniqueness of the equilibrium point, we only 
need to prove that (|9|l has a unique solution. 

To this end, we check equation ^ element-wisely. Making 
use of the structure of matrix A in (|6]l, we have the following 
equation for alH = 1, 2, . . . , A^; 

(10) 
Since the interaction is bi-directional, i.e., a^.j = aj,i, it 
follows 

N N 

Next, we prove that (fT2] i holds by proving that both 
max{g(x3 J} and min{_g(x3 J} are zero; 



gi^li) = 5(2^3.2) = • • ■ = 9{xIn) ^ 



(12) 



Suppose to the contrary that (fT2t does not hold, and hence 

nia.x{g{x^ J} > is satisfied since J2i=i 5(^3 i) = holds 
according to ([Til l. Represent the index that has the largest 
g{xl j) among all z (i = 1, 2, . . . , N) as m. If there are mul- 
tiple indices m satisfying g{xl „j) = max{gr(a:3 j)}, then any 

' i ' 

one can be index m. Because /(•) is a decreasing function, 
it follows that g{») is a decreasing function according to its 
definition on the left hand side of fTOl i. So if the maximal 
.9(^3 i) i^ attained when 



among 2:5 1, XI2, 
of ([Toll, i.e.. 



should be the smallest 
,2:3^. Therefore, the right hand side 



■'3,m 



bib3y^amj{x: 






(13) 



should be non-positive, and hence 5(2:3^) should be non- 
positive. This contradicts the fact that g{xl „) is the largest 
among all g{xl J and it is positive (due to the constraint in 
(ITTli). Hence ina.x{g{xl J} = holds. 

i ■ 

Similarly, we can prove that min{5(a;3 J} — holds. 

i 

Therefore, we have il2t . which further leads to 



/(a^3,,) = bib2b3xl.^, i = l,2,. 



,N 



(14) 



Given that /(•) is a monotonic decreasing function on M+, 
it follows that the solution to ([T4l i is unique and it satisfies 

■^3 1 =X*32 = ■■■ = -^3 N ^ ^0 > 0, f{xo) = 616263X0 

(15) 

Therefore the solution to (|9]l is unique, thus the equilibrium 
point is unique and hence isolated. 

Step II — Satisfaction of condition (b): Following the 
derivations in [24], [25], we can easily get that condition 
(b) is satisfied. 



Step in — Satisfaction of condition (c): The fact that J has 
at least one eigenvalue with positive real part is equivalent 
to the statement that the linearized system of dU around the 
equilibrium point is strictly unstable. So instead of proving 
(c) directly, next we prove the strict instability of linearized 
system of (01) around the equilibrium point under condition 
(|7]i. To this point, we transform (|4) into the frequency domain 
as shown in Fig. |2l where H(s) is given by 

H{s) = {{si + hil){sl + 62/ + A){sl + 63/))"^ 

1 1 (16) 

{s + bi){s + h^) 



His) 



Fig. 2. Schematic diagram of tlie frequency domain formulation of 0. 

Linearizing the nonlinear item /(X3) in (|4|i around the 

xq ] yields 



equilibrium point X^ 



[ Xq Xn 

{i+4r 



PXq {bib2b3Xoy 



(17) 

In the second equality, the relation f{xo) = bib2b3Xo at the 
equilibrium point is employed. 

Based on the linearization in (fTTI i. we can get the closed- 
loop dynamics of the system in Fig. [2] as 



G{s) = {I-<jH{s))-'H{s) 



(18) 



Since ^ is a symmetric matrix, it only has real eigenvalues 
and it can always be diagonalized as follows: 

A = PTP-\ T = diag(ui,^;2,.-.^^iv) (19) 

where P is the similarity transformation matrix. Since A is 
the graph Laplacian, it always has an eigenvalue associated 
with an eigenvector composed of identical elements [26]. 
Here we arrange the eigenvalues in increasing order, so 
ui = always holds and the connectivity assumption in 
Assumption [U leads to Vi > (i = 2,S, . . . , N) [26]. 
Substituting ( fT9l ) into ( fT6l ). we have 

H{.s) = PKP-\ A = diag(Ai, A2, . . . A^) (20) 

where eigenvalues A; (i = 1, 2, . . . , N) are given by 

Ai = 



1 



, j=2,3,...,iV 



(s + &i)(s + 52)(.s + 63) 

^ ^ (s + fel)(s + 62+Wj)(s + 63) 

Using ( fTsT i and ( |20] |. we can obtain that 

G(s) = P(/ - crA)-iAP-i = PAP-i (21) 



where 

A = diag((5i, ^2,..., 5n) 

with eigenvalues 5i (i = 1, 2, . . . , N) given by 
Ai 1 



5i = 



1-crAi {s + bi){s + b2){s + b:^) - a 



<S, = 



1 - aXj {s + bi){s + 62 + V3){s + 63) - cr' 



(22) 



(23) 



where j = 2, 3, . . . , A^. 

According to the Routh— Hurwitz stability criterion, we 
know 5i is strictly unstable if and only if 

a < 616263 " (61 + 62 + 63)(6i62 + 6163 + 6263) (24) 



is satisfied, and Sj {j = 2,3,. 
and only if 



. , A^) is strictly unstable if 



(25) 



a <6i(62 + ■L'j)63 - (61 + 62 + Vj + 63) X 
(61(62 + Vj) + 6163 + (62 + Wj)63) 

is satisfied. Note that JTM is a necessary condition for 
(I25I 1. so G{s) is strictly unstable if and only if (|24] | holds. 
Substituting a in ( [TtI i into (|24] |. we know that G{s) is strictly 
unstable if and only if the inequality (|7]) in Theorem [1] holds, 
i.e., condition (c) holds if (jT) is satisfied. 

So far, we have proven conditions (a), (b), and (c), and 
hence guaranteed the existence of oscillatory solutions to (|4]l. 
The synchronization condition for these oscillations, i.e., (jSj, 
follows easily from the secant condition in [14], thus the 
proof is omitted due to space limitations. Hence Theorem [T] 
is proven. ■ 

Remark 3: From ^, we can see that with an increase in 
7, a stronger network connectivity g is needed to guarantee 
network synchronization. Given that for p > 1, j can be 
verified an increasing function of the Hill coefficient p, 
we know that a system having a higher Hill coefficient 
(i.e., higher end product cooperativity) requires a stronger 
coupling to maintain network synchronization. 

IV. Oscillation period estimation based on 

MULTIVARIABLE HARMONIC BALANCE TECHNIQUE 

A. Oscillation analysis based on harmonic balance tech- 
nique 

In this section, we reformulate the problem of oscillation 
analysis using a multivariable harmonic balance technique 
[20]. This is motivated by the observation that II{s) is a 
low pass filter thus higher order harmonics of oscillations 
in the close-loop system can be safely neglected. Hence the 
waveform of 2:3 ^ (i — 1,2, ... ,N) can be approximated 
by its zero-order and first-order harmonic components [20], 
[27]: 



X3A = at + Pi s\\i{wt + 4>i), i = 1,2,. . .,N 



(26) 



where ai and /3i denote the amplitudes of the zero-order and 
the first-order harmonic components, respectively, and w and 
(j)i denote the oscillation frequency and phase, respectively. 



Since /(•) is a static nonlinear function, it can be approx- 
imated by its describing functions [27]: 



where 



/(a;3,j) ~ itai + ViPi sin(wi + (i>{) 



f{a, + 13, sin{t))dt 



27ra,; 



Vi = 



27ra,- 



/(ofi + Pi sin(i)) sm{t)dt 



(27) 



(28) 



(29) 



The describing function ^j is the gain of /(•) when the input 
is a constant value a; and the output is approximated by its 
zero-order harmonic. The describing function rji is the gain 
of /(•) when the input is a sinusoid of amplitude /3i and the 
output is approximated by its first-order harmonic [27]. 

Consequently, the closed-loop equations that ai and f3i are 
expected to satisfy are given by [20] 



{I -H(Q)E)a = 



and 



(/ - H{jw)n)l3 = 



(30) 



(31) 



respectively, where 



ai 



,13 = 



PieJ')'^ 



13 NC- 



J<PN 



^ = diag{^i, ...,^n}, 
n = diag{77i, ...,f]N} 



Therefore, the problem of oscillation analysis reduces to 
finding a, (3, and w satisfying dSOl l and dSTT ). Note that ( l30b 
and dSTI ) are referred to as harmonic balance equations. 

Let S* and H* be constant matrices satisfying dSOl l and 
dJTT i simultaneously. Define two Unear systems Go(s) and 
6*1(5) as 



Giis)^iI-H{s)U*)-'His) 



(32) 



The systems Gois) and Gi(s) are obtained by replacing 
the nonlinearity /(•) with the constant gain computed from 
the describing functions. Thus, the two linear systems con- 
tain some information about the oscillations of the original 
nonlinear system. According to Iwasaki [20], the predicted 
oscillation at frequency oj is expected stable if both Go(s) 
and Gi(s) are marginally stable with poles of s = and 
s = ijw on the imaginary axis, respectively (the rest in the 
open left half plane). Therefore, the problem of oscillation 
analysis can be reduced to the following problem: 

Problem 1: For the coupled Goodwin oscillators in (|4]i, 
find S* and 11* that 

• satisfy ( l30l i and ( 1311 ). respectively; 

« and at the same time guarantee that Go(s) and Gi(s) 

in ( |32] | are marginally stable. 
The solution is given in the next section. 



B. Oscillation period of coupled Goodwin oscillators 

According to [20], ( l30l l and dSTT i are very difficult to solve 
since in general S and 11 depend on a and j3. Bearing in 
mind that we are interested in the collective period, we can 
restrict our attention to solutions that describe synchronized 
oscillations of the oscillator network. This provides a clue 
to solve the problem: according to Definition |2] synchrony 
means that x-^^i are identical, i.e., 1) the phases (pi (i = 
1,2,..., N) are identical; 2) the amplitudes a^ and (3i {i = 
1,2,..., N) are identical. Given that ^i and rji are determined 
by tti and /3i, we further have the equality of all ^^ and all 
rji. Making use of these properties, we have Theorem |2] 

Theorem 2: For the Goodwin oscillator network in (2), 
if the oscillation/synchronization condition in Theorem [T] is 
satisfied, then its collective period Tcoiicctivc is given by 

27r 27r 



Tr. 



)llccti' 



Proof: According to the above analysis, we have 

d = al, (3 — /31, S — ^In, n = ry/jv 
where a, f3, ^, and 77 are constants and 1 is given by: 



(33) 



(34) 



1=[1 1 
Hence ( [30b and dST] ) reduce to 



1]^ 



and 



(-I-H{0))al 



(-1 - H(iw))(3l = 



(35) 



(36) 



respectively, which further means that j is the eigenvalue 
of H{Q) corresponding to the eigenvector with identical 
elements, and - is the eigenvalue of H{jw) corresponding 
to the eigenvector with identical elements. 

From (|20] |. we know the eigenvalues of H{Q) are Ai = 
, } , and A," = ,. ,, ] — rrr for 7 = 2, 3, . . . , A^. Further 
notice that only Ai corresponds to eigenvectors with identical 
elements. Thus we have 



\- 



1 



hhh 



(37) 



Similarly, we can get that the eigenvalues of 



Hijw) 



are 



Al 



{jw+bi) {jw+b2 ){jw+b3) 

for J = 2,3,..., N. 



)^ . _ i 

1 {j'W+bi ) (jw+b2 +Vj ) ( JMi+63 ) 

Further notice that only Ai corresponds to eigenvectors with 
identical elements. Thus we have 



1] {jw + bi){jw + b2)ijw - 



(38) 



According to (|29l l, r/ is real, thus the item on the right hand 
side of ( |38] | must be real, i.e., its imaginary part is zero. 
Given that 

{jw + bi){jw + b2){jw + 63) = 61&2&3 - w'^{bi + 62 + 63) 

+ jw((bib2 + bibs + b2b3) - w^) 

(39) 



we have the imaginary part jw{{bib2 + foi&a + fo2&3) — w^) 
equal to 0, which further leads to 

,.2 



W 



6l&2 + ^1^3 +62^3 



(40) 



and 



■q = 61&2&3 - (&1&2 + &1&3 + &2&3)(&1 + &2 + ^3) (41) 

Hence we know the collective period is determined by (l33T l. 

But to prove Theorem |2] it remains to prove that the 
oscillation at estimated frequency is stable, i.e., Gq{s) and 
Gi(s) in (|32] i are marginally stable [20]. So we proceed to 
prove that: (1) Go{s) has one pole of s = and the rest in the 
open left half plane and, (2) Gi (s) has a pair of imaginary 
poles s — ±jw and the rest in the open left half plane. 

Substituting S and 11 in ( [34l i into ( [32l ) yields 



Go(,s) = (/-eii^(,s))-iH(s) 



and 



(42) 
(43) 



with ^ and 77 given in (|37] | and (|4TJ, respectively. 

We first consider Go(s) in (gill. From ^, ^, and (|23]l 
we know the eigenvalues of Go(s) are given by 

1 



Si = 



(S + bl)is + b2)is + bs) - e 

1 



"' {s + b,)is + b2 + v,){s + bs)-C ^ 2,3,. ..,Ar 

Substituting ^ in dJTJ i into the above equations, we know the 
poles of ( |42] | are determined by the roots of 

(s + 6i)(s + 62)(s + &3)-&ifc263 = (44) 

and 



(s + 6i)(s + &2 + t^i)(s + 63) - 616263 = 



(45) 



forj=2,3,...,iV. 

It is clear that (|44] | has one root s = 0. And using the 
Routh— Hurwitz stability criterion, we can get that the other 
roots of (|44] | and all roots of (|45t are in the open left half 
plane. Hence Go(s) is marginally stable. 

Similarly, we can prove that the eigenvalues of Gi (s) are 



and 



5, 



(s + 6i)(s + 62)(s + 63)-r? 
1 



(.s + 6i)(s + 62+i;j)(s + 63)-r; 

for j = 2,3,...,A^, with -q given in ( 1411 1. And hence its 
poles are determined by the roots of 

(s + 6i)(s + 62)(s + 63) - 616263 

+ (6162 + 6163 + 6263)(6i + 62 + 63) 

and 

(s + 6i)(s + 62 + Wj)(s + 63) - 616263 

+ (6162 + 6163 + 6263)(6i + 62 + 63) 

forj=2,3,...,iV. 











(46) 



(47) 



It can be verified that s = ±jw are two roots of 
And using the Routh— Hurwitz stability criterion, it follows 
that the other roots of (|46] | and all roots of (l47t are in the 
open left half plane. So Gi(s) is marginally stable. Hence 
the derived oscillation at frequency w in ( |33] ) is stable, which 
completes the proof. ■ 

Remark 4: The collective period in ( l33T l is given in terms 
of the dimensionless parameters in (l2). Representing the 
collective period with the original dimensional parameters 

gives 

_ 27r 

-^ collective — / , , , . , ^= \^'~>) 

\Jk\k2 + fclK3 + K.2k-i 

Eqn. ( |48] l means that the collective period decreases with 
an increase in the rate constants for degradation of each 
component, but it is independent of the rates of transcription, 
translation, and catalysis. These give insights into the basic 
determination mechanism of the collective period in coupled 
biological oscillator networks, and may further provide guid- 
ance in synthetic biology design. 

Remark 5: From (l48l l. we can see that when the intercel- 
lular interaction is of the form in (|3), the collective period 
is only determined by k\, fc2, and ^3, and it is independent 
of intercellular coupling. The results are obtained based on 
analytical treatment of a network of coupled gene regulatory 
oscillators and they corroborate the results in [19], which are 
obtained using the phenomenological single-variable phase 
model and state that the strength of intercellular coupling 
does not affect the collective period of circadian rhythm 
oscillator networks. 

Remark 6: It is worth noting that if the coupling is of 
a form different from (|3), then it may affect the collective 
period, even if it is still of diffusive type. Examples have 
been reported in [2] and [3]. 

V. Numerical study 

In this section, simulation results are given to confirm 
the theoretical predictions. We considered a network of nine 
Goodwin oscillators. The coupling strengths a^ j were ran- 
domly chosen from a uniform distribution on [0, 1] and are 
given in Table I. It can be verified that the coupling topology 
is connected and the algebraic connectivity \s g ~ 2.4583 
({3 is equal to the second smallest eigenvalue of interaction 
matrix A in ^, as defined in Theorem [T]i. 

TABLE I 
Coupling TOPOLOGY ai_j (1 < «,i < 9, « ^ j) of the Goodwin 

OSCILLATOR NETWORK 



^\ 


1 


2 


3 


4 


5 


6 


7 


8 


9 


1 




0.3 


0.5 





0.6 


0.2 





0.7 


0.8 


2 


0.3 




0.7 


0.2 


0.1 


0.8 


0.3 


0.1 


0.5 


3 


0.5 


0.7 




0.3 


0.6 


0.2 


0.6 





0.8 


4 





0.2 


0.3 




0.4 


0.6 


0.2 


0.9 


0.1 


5 


0.6 


0.1 


0.6 


0.4 




0.2 


0.7 


0.3 


0.8 


6 


0.2 


0.8 


0.2 


0.6 


0.2 




0.1 


0.9 


0.3 


7 





0.3 


0.6 


0.2 


0.7 


0.1 




0.4 


0.5 


8 


0.7 


0.1 





0.9 


0.3 


0.9 


0.4 




0.8 


9 


0.8 


0.5 


0.8 


0.1 


0.8 


0.3 


0.5 


0.8 





First we tested our oscillation/synchronization condition 
in Theorem [T] The resuhs are summarized in Table II. It can 
be seen that oscillation/synchronization can be obtained only 
when the parameters satisfy i? > 1 in (|7]). 

TABLE II 

Test of the oscillation/synchronization condition 



p 


bi 


b2 


bz 


R 


Simulation results 


17 


0.4 


0.4 


0.4 


1.7102 


Oscillation/synchronization 


17 


0.5 


0.5 


0.5 


1.6541 


Oscillation/synchronization 


17 


0.6 


0.6 


0.6 


1.5286 


Oscillation/synchronization 


17 


0.7 


0.7 


0.7 


1.3266 


Oscillation/synchronization 


17 


0.8 


0.8 


0.8 


1.0421 


Oscillation/synchronization 


17 


0.85 


0.85 


0.85 


0.8686 


No oscillation/synchronization 


17 


0.9 


0.9 


0.9 


0.676 


No oscillation/synchronization 


17 


1.0 


1.0 


1.0 


0.2620 


No oscillation/synchronization 


17 


0.7 


0.8 


0.9 


1.0433 


Oscillation/synchronization 


17 


0.9 


0.8 


0.8 


0.9300 


No oscillation/synchronization 



We also compared the collective periods in the oscillatory 
cases in Table II with the estimated value. The results are 
summarized in Table III. It can be seen that the estimated 
values approximate the actual collective periods very closely. 

TABLE III 

Comparison between the estimated collective period [s] and 

the actual collective period [s] 



61=62= 63 


0.4 


0.5 


0.6 


0.7 


0.8 


Actual value 


10.68 


8.00 


6.31 


5.23 


4.53 


Estimated value 


11.35 


7.26 


6.05 


5.19 


4.54 


Estimation error 


6.27% 


-9.25% 


-4.12% 


-0.76%- 


0.22% 



VI. Conclusions 

Underlying biological rhythms are networks of interacting 
cellular oscillators. How the collective oscillation patterns 
arise from autonomous cellular oscillations is poorly under- 
stood. The Goodwin oscillator is a quintessential model of 
biochemical oscillators based on negative feedback mecha- 
nisms. Based on a network of coupled Goodwin oscillators, 
we studied analytically the oscillation/synchronization con- 
dition and collective period of coupled biochemical oscilla- 
tors by using a multivariable harmonic balance technique. 
We give an oscillation/synchronization condition of coupled 
Goodwin oscillators. The condition shows that a system 
having a higher Hill coefficient (corresponding to a higher 
cooperativity of end product repression) requires a stronger 
intercellular coupling to maintain network synchronization. 
We also analytically estimate the collective oscillation period 
of the oscillator network. The collective oscillation period is 
only dependent on the degradation rates of each component, 
and it is independent of the rate of transcription, translation, 
and catalysis. The results are confirmed by numerical sim- 
ulations and may provide guidance in synthetic-biological- 
oscillator design. Given that the Goodwin oscillator has 
been successfully implemented in vivo, synthetic biology 
based testing of the predictions is promising. Experimental 
verification is also feasible in biological oscillators whose 
degradation/synthesis rates are tunable. 
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